function YY = CGH_steadystate_solve(x)

global M_ normalizationI fixed_costI capitalI

% read out parameters to access them with their name
NumberOfParameters = M_.param_nbr;
for ii = 1:NumberOfParameters
  paramname = deblank(M_.param_names(ii,:));
  eval([ paramname ' = M_.params(' int2str(ii) ');']);
end


Y       = x(1);
alppha  = x(2);
eta     = x(3);
alpha1P = x(4); 
alpha2P = x(5); 

load indicators

N_n              = 0.19*homeprodI + (1-homeprodI)*1e-10; 
N_m              = 0.33; 
K_m              = 5.15*Y*capitalI + 1 - capitalI;
K_n              = 6.75*Y*home_capitalI*capitalI + 1 - home_capitalI;

M                = betta;
mu               = theta_mu/(theta_mu-1);
Xi               = 1/(subsidy*mu);
Z                = 1;
u                = 1;
u_n              = 1;
q_n              = 1;
q                = 1;
delta_1          = 1/betta-1+delta_0;
delta_n_1        = 1/betta-1+delta_n_0;
R_R              = 1/betta;

if capitalI==1
    PF_normalization = normalizationI*(1/(Xi*N_m^(1-alppha)*K_m^alppha)) + (1-normalizationI);
    Psi              = fixed_costI*(1-Xi)*PF_normalization*N_m^(1-alppha)*K_m^alppha;
    I_m              = delta_0*K_m;
    I_n              = delta_n_0*K_n*home_capitalI;
    W                = ((1-alppha)*Xi*PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha))/N_m;
    R_K              = (alppha*Xi*PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha))/(u*K_m);
    D                = Y - N_m*W - delta_0*K_m;   
    B                = nu*K_m;
    D_E              = D-(1-1/R_R)*B;
    P_E              = betta/(1-betta)*D_E;
    C_m              = - I_n + - (P_E+1/R_R*nu*K_m) + W*N_m + D_E + P_E + nu*K_m;
else
    PF_normalization = normalizationI*(1/(Xi*N_m)) + (1-normalizationI);
    Psi              = fixed_costI*(1-Xi)*PF_normalization*N_m;
    I_m              = 0;
    I_n              = 0;
    W                = Xi*PF_normalization*(Z*N_m)/N_m;
    D                = Y - N_m*W;   
    R_K              = 1;
    B                = 0;
    D_E              = D;
    P_E              = 0;
    C_m              = - I_n + W*N_m;
end

if homeprodI==1
    if home_capitalI==1
       C_n              = ((u_n*K_n)^alpha2P)*(N_n^(1-alpha2P));
    else
       C_n              = N_n;    
    end
    C                = (alpha1P*(C_m^b1P) + (1-alpha1P)*(C_n^b1P))^(1/b1P) ;
else
    C_n = 1e-10;
    C                = C_m ;
end

L                = 1 - N_m - N_n;
utilC            = eta*(C^(eta-1))*(L^(1-eta));
utilL            = (1-eta)*(C^eta)*(L^-eta);


if homeprodI==1
   if home_capitalI==1 && capitalI==1
      YY = [ - Y + PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha)-Psi;
             - W + (utilL/utilC)*(1/alpha1P)*((C_m/C)^(1-b1P));      
             - q + M*(u*R_K+q*(1-(delta_0 + delta_1*(u-1)+delta_2/2*(u-1)^2)-phi_k/2*(I_m/K_m-delta_0)^2 +phi_k*(I_m/K_m-delta_0)*(I_m/K_m)));
             - utilL/utilC + (1-alpha1P)*(1-alpha2P)*((C_n/C)^(b1P - 1))*(C_n/N_n); 
             - q_n*u_n*delta_n_1 + alpha2P*((1-alpha1P)/alpha1P)*(C_n/K_n)*((C_n/C_m)^(b1P-1));];  
   elseif capitalI==0
      YY = [ - Y + PF_normalization*(Z*N_m)-Psi;
             - W + (utilL/utilC)*(1/alpha1P)*((C_m/C)^(1-b1P));      
             - alppha + 1e-12;
             - utilL/utilC + (1-alpha1P)*(1-alpha2P)*((C_n/C)^(b1P - 1))*(C_n/N_n); 
             - alpha2P;];  
   else
     YY = [ - Y + PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha)-Psi;
            - W + (utilL/utilC)*(1/alpha1P)*((C_m/C)^(1-b1P));      
            - q + M*(u*R_K+q*(1-(delta_0 + delta_1*(u-1)+delta_2/2*(u-1)^2)-phi_k/2*(I_m/K_m-delta_0)^2 +phi_k*(I_m/K_m-delta_0)*(I_m/K_m)));
            - utilL/utilC + (1-alpha1P)*((C_n/C)^(b1P - 1))*(C_n/N_n);
            - alpha2P];
   end
else
    if capitalI==1
        YY = [ - Y + PF_normalization*(u*K_m)^alppha*(Z*N_m)^(1-alppha)-Psi;
               - W + (utilL/utilC)*(1/alpha1P)*((C_m/C)^(1-b1P));      
               - q + M*(u*R_K+q*(1-(delta_0 + delta_1*(u-1)+delta_2/2*(u-1)^2)-phi_k/2*(I_m/K_m-delta_0)^2 +phi_k*(I_m/K_m-delta_0)*(I_m/K_m)));
               - alpha1P + 1;
               - alpha2P]; 
    else
        YY = [ - Y + PF_normalization*(Z*N_m)-Psi;
               - W + (utilL/utilC)*(1/alpha1P)*((C_m/C)^(1-b1P));      
               - alppha + 1;
               - alpha1P + 1;
               - alpha2P]; 
    end
    
end